
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Test: Choice of Specification for the Bidders' Private Values    %
%                      Linear versus Exponential                         %
%                         Equivalence Property                           %
%                        Monte Carlo Experiment                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all

L2=200;
L3=200;
L=L2+L3;
N=[2,3];
alpha=0.50:0.10:0.80; % quantile levels alpha
psi2=ones(L2,1)*(N(1)*(alpha.^(N(1)-1))-(N(1)-1)*(alpha.^N(1))); % Psi function given N=2
psi3=ones(L3,1)*(N(2)*(alpha.^(N(2)-1))-(N(2)-1)*(alpha.^N(2))); % Psi function given N=3
psi = [psi2;psi3];
m=size(alpha,2);

eps2=unifrnd(0,1,L2*N(1),1); % Regression Errors N=2
eps2=reshape(eps2,N(1),L2);
eps2=sort(eps2,'descend');
eps2=eps2(2,:);
eps3=unifrnd(0,1,L3*N(2),1); % Regression Errors N=3
eps3=reshape(eps3,N(2),L3);
eps3=sort(eps3,'descend');
eps3=eps3(2,:);
eps=[eps2';eps3'];

x2=[ones(L2,1),normrnd(0,1,L2,1)]; % Covariates N=2
x3=[ones(L3,1),normrnd(0,1,L3,1)]; % Covariates N=3
x=[x2;x3];
p=size(x,2);

gammai=[0,0.5]'; % Initial guess for coefficients 

YL=x*gammai+eps; % Winning bid specification under a linear model for private values
YE=exp(x*gammai+eps); % Winning bid specification under a exponential model for private values

% Auxiliary parameters
toler=0.1;
iter=1000;
Lambda=linspace(0,1,iter);
Delta=zeros(iter,m);

% Search of lambda (in Lambda) that satisfies equivalence property between Linear and exponential specifications
for i=1:iter;
  lambda=Lambda(i);  
  y=lambda*YL+(1-lambda)*YE; % linear combination
    for k=1:m
       [objL ~]=MM(psi(:,k), p, y, x, toler, gammai); % objective function of the linear model
       [objE ~]=MM_Exp(psi(:,k), p, y, x, toler, gammai); % objective function of the exponential model
       Deltak(:,k)=(objE-objL)/objL; % relative difference between both optimized objective functions
    end

  Delta(i,:)=Deltak;
  i
end


for k=1:m
    Deltamin = min(abs(Delta(:,k)));
    [~,indx]=ismember(Deltamin,abs(Delta(:,k)));
    lambdastar(:,k)=Lambda(:,indx); % lambda that minimizes Deltak for each k
    Dmin(:,k)=Delta(indx,k);
end

Deltasup=max(Dmin)
